clc;clear;
tic;
% hetero-nuclear PHIP program for 6 spins
format long e

% evaluating the energy levels of the 6 spin system
ai=sqrt(-1);

sx=[0 0.5; 0.5 0];                  % matrix of S_x for spin 1/2
sy=[0 -0.5*ai; 0.5*ai 0];           % matrix of S_y for spin 1/2
sz=[0.5 0; 0 -0.5];                 % matrix of S_z for spin 1/2
od=[1 0; 0 1];                      % 2*2 unity matrix

% evaluating the energy levels of the 6 spin system
nspins=6;
matrixsize=2^6;

% setting {x,y,z} components of spin operators of all spins: 1024*1024
% matrices
I1z=zeros(64,64);
I2z=zeros(64,64);
I3z=zeros(64,64);
I4z=zeros(64,64);
I5z=zeros(64,64);
I6z=zeros(64,64);

I1x=zeros(64,64);
I2x=zeros(64,64);
I3x=zeros(64,64);
I4x=zeros(64,64);
I5x=zeros(64,64);
I6x=zeros(64,64);


I1y=zeros(64,64);
I2y=zeros(64,64);
I3y=zeros(64,64);
I4y=zeros(64,64);
I5y=zeros(64,64);
I6y=zeros(64,64);


% setting matrices of scalar products of all spins of interest
pro12=zeros(64,64);
pro14=zeros(64,64);

pro32=zeros(64,64);
pro34=zeros(64,64);

pro52=zeros(64,64);
pro54=zeros(64,64);

pro26=zeros(64,64);
pro46=zeros(64,64);


% pro13=zeros(1024,1024);
% pro14=zeros(1024,1024);
% pro15=zeros(1024,1024);
% pro23=zeros(1024,1024);
% pro24=zeros(1024,1024);
% pro25=zeros(1024,1024);
% 
% pro16=zeros(1024,1024);
% pro17=zeros(1024,1024);
% pro26=zeros(1024,1024);
% pro27=zeros(1024,1024);
% pro67=zeros(1024,1024);
% 
% pro68=zeros(1024,1024);
% pro79=zeros(1024,1024);
% pro69=zeros(1024,1024);
% pro78=zeros(1024,1024);
% 
% pro610=zeros(1024,1024);
% pro710=zeros(1024,1024);
% pro810=zeros(1024,1024);
% pro910=zeros(1024,1024);


Ix=zeros(64,64);
Iy=zeros(64,64);
Iz=zeros(64,64);


I1z=kron(sz,kron(od,kron(od,kron(od,kron(od,od)))));
I2z=kron(od,kron(sz,kron(od,kron(od,kron(od,od)))));
I3z=kron(od,kron(od,kron(sz,kron(od,kron(od,od)))));
I4z=kron(od,kron(od,kron(od,kron(sz,kron(od,od)))));
I5z=kron(od,kron(od,kron(od,kron(od,kron(sz,od)))));
I6z=kron(od,kron(od,kron(od,kron(od,kron(od,sz)))));

I1x=kron(sx,kron(od,kron(od,kron(od,kron(od,od)))));
I2x=kron(od,kron(sx,kron(od,kron(od,kron(od,od)))));
I3x=kron(od,kron(od,kron(sx,kron(od,kron(od,od)))));
I4x=kron(od,kron(od,kron(od,kron(sx,kron(od,od)))));
I5x=kron(od,kron(od,kron(od,kron(od,kron(sx,od)))));
I6x=kron(od,kron(od,kron(od,kron(od,kron(od,sx)))));

I1y=kron(sy,kron(od,kron(od,kron(od,kron(od,od)))));
I2y=kron(od,kron(sy,kron(od,kron(od,kron(od,od)))));
I3y=kron(od,kron(od,kron(sy,kron(od,kron(od,od)))));
I4y=kron(od,kron(od,kron(od,kron(sy,kron(od,od)))));
I5y=kron(od,kron(od,kron(od,kron(od,kron(sy,od)))));
I6y=kron(od,kron(od,kron(od,kron(od,kron(od,sy)))));

% 
% I1z=kron(sz,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I2z=kron(od,kron(sz,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I3z=kron(od,kron(od,kron(sz,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I4z=kron(od,kron(od,kron(od,kron(sz,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I5z=kron(od,kron(od,kron(od,kron(od,kron(sz,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I6z=kron(od,kron(od,kron(od,kron(od,kron(od,kron(sz,kron(od,kron(od,kron(od,od)))))))));
% I7z=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sz,kron(od,kron(od,od)))))))));
% I8z=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sz,kron(od,od)))))))));
% I9z=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sz,od)))))))));
% I10z=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,sz)))))))));
% 
% I1x=kron(sx,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I2x=kron(od,kron(sx,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I3x=kron(od,kron(od,kron(sx,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I4x=kron(od,kron(od,kron(od,kron(sx,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I5x=kron(od,kron(od,kron(od,kron(od,kron(sx,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I6x=kron(od,kron(od,kron(od,kron(od,kron(od,kron(sx,kron(od,kron(od,kron(od,od)))))))));
% I7x=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sx,kron(od,kron(od,od)))))))));
% I8x=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sx,kron(od,od)))))))));
% I9x=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sx,od)))))))));
% I10x=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,sx)))))))));
% 
% I1y=kron(sy,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I2y=kron(od,kron(sy,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I3y=kron(od,kron(od,kron(sy,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I4y=kron(od,kron(od,kron(od,kron(sy,kron(od,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I5y=kron(od,kron(od,kron(od,kron(od,kron(sy,kron(od,kron(od,kron(od,kron(od,od)))))))));
% I6y=kron(od,kron(od,kron(od,kron(od,kron(od,kron(sy,kron(od,kron(od,kron(od,od)))))))));
% I7y=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sy,kron(od,kron(od,od)))))))));
% I8y=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sy,kron(od,od)))))))));
% I9y=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(sy,od)))))))));
% I10y=kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,sy)))))))));


p12=I1x*I2x+I1y*I2y+I1z*I2z;
p14=I1x*I4x+I1y*I4y+I1z*I4z;

p32=I3x*I2x+I3y*I2y+I3z*I2z;
p34=I3x*I4x+I3y*I4y+I3z*I4z;

p52=I5x*I2x+I5y*I2y+I5z*I2z;
p54=I5x*I4x+I5y*I4y+I5z*I4z;

p26=I2x*I6x+I2y*I6y+I2z*I6z;
p46=I4x*I6x+I4y*I6y+I4z*I6z;

% p13=I1x*I3x+I1y*I3y+I1z*I3z;
% p14=I1x*I4x+I1y*I4y+I1z*I4z;
% p15=I1x*I5x+I1y*I5y+I1z*I5z;
% p23=I2x*I3x+I2y*I3y+I2z*I3z;
% p24=I2x*I4x+I2y*I4y+I2z*I4z;
% p25=I2x*I5x+I2y*I5y+I2z*I5z;
% 
% p16=I1x*I6x+I1y*I6y+I1z*I6z;
% p17=I1x*I7x+I1y*I7y+I1z*I7z;
% p26=I2x*I6x+I2y*I6y+I2z*I6z;
% p27=I2x*I7x+I2y*I7y+I2z*I7z;
% p67=I6x*I7x+I6y*I7y+I6z*I7z;
% 
% p68=I6x*I8x+I6y*I8y+I6z*I8z;
% p79=I7x*I9x+I7y*I9y+I7z*I9z;
% p69=I6x*I9x+I6y*I9y+I6z*I9z;
% p78=I7x*I8x+I7y*I8y+I7z*I8z;
% 
% p610=I6x*I10x+I6y*I10y+I6z*I10z;
% p710=I7x*I10x+I7y*I10y+I7z*I10z;
% p810=I8x*I10x+I8y*I10y+I8z*I10z;
% p910=I9x*I10x+I9y*I10y+I9z*I10z;
% 
% p18=I1x*I8x+I1y*I8y+I1z*I8z;
% p28=I2x*I8x+I2y*I8y+I2z*I8z;
% p19=I1x*I9x+I1y*I9y+I1z*I9z;
% p29=I2x*I9x+I2y*I9y+I2z*I9z;
% 
% p110=I1x*I10x+I1y*I10y+I1z*I10z;
% p210=I2x*I10x+I2y*I10y+I2z*I10z;



Iz=I1z+I2z+I3z+I4z+I5z+I6z;
Ix=I1x+I2x+I3x+I4x+I5x+I6x;
Iy=I1y+I2y+I3y+I4y+I5y+I6y;



% J-couplings in Hz
J12=7.5;
J14=7.5;

J32=7.5;
J34=7.5;

J52=7.5;
J54=7.5;

J26=1.1;
J46=1.1;

% J13=7.55;
% J14=7.55;
% J15=7.55;
% J23=7.55;
% J24=7.55;
% J25=7.55;
% 
% J16=-0.6;
% J17=-0.6;
% J26=-0.6;
% J27=-0.6;
% J67=1.5;
% 
% J68=7.5;
% J79=7.5;
% J69=0.8;
% J78=0.8;
% 
% J610=1.3;
% J710=1.3;
% J810=7.5;
% J910=7.5;
% 
% J18=0.29;
% J28=0.29;
% J19=0.29;
% J29=0.29;
% 
% J110=-0.45;
% J210=-0.45;




tp=10.0;   % preparation time of PHIP: should be long that all coherences are washed out
te=10.0;     % free evolution time of PHIP at the preparation field, should be usually set to zero

% all chemical shifts in ppm in the sigma-scale
ch1=-0.95;
ch3=-0.95;
ch5=-0.95;

ch2=-1.65;
ch4=-1.65;

ch6=-9.44;


% Magnetic fields
% FV: time and number of points:
T2=2.0; % T_2 value, which gives the NMR linewidth 

[BotT] = textread('Bott01mT.DAT');  % B(t) profile of field switching
Nfv=2995;                           % number of points in the profile
        
B=10000*BotT(1,2);                  % starting (preparation) field
Bdet=70400;                         % final (detection) field

H = zeros(64,64);               % Hamiltonian: 64/64 matrix
    rho = zeros(64,64);         
    rho0 = zeros(64,64);
    rhot = zeros(64,64);
    
    rhos=zeros(4,4);
    rhos(2,2)=1/2;
    rhos(3,3)=1/2;
    rhos(2,3)=-1/2;
    rhos(3,2)=-1/2;
    

    H=B*3.0013e2*(ch1*I1z+ch2*I2z+ch3*I3z+ch4*I4z+ch5*I5z+ch6*I6z)/7.04e4;
    H=H+J12*p12+J14*p14+J32*p32+J34*p34+J52*p52+J54*p54+J26*p26+J46*p46;
    
% setting starting density matrix: 64*64; spins 1 and 2 are ALWAYS
% coming from para-hydrogen

rho0=kron(rhos,kron(od,kron(od,kron(od,od))));
rho0=rho0/(2^(nspins-2));

% rho0=kron(od,kron(rhos,kron(od,kron(od,kron(od,kron(od,kron(od,kron(od,od))))))));
%         rho0=rho0/256;
        

% calculating the transition matrix 


    

% calculating the eigenvenctors
evec=zeros(64,64);
evec1=zeros(64,64);
energ=zeros(64,64);
HH=zeros(64,64);
[evec,energ]=eig(H);    % solving the eigen-problem for the Hamiltonian




evec;
evec1=inv(evec);  
        HH=evec1*H*evec;
        
        rho0=evec1*rho0*evec;   % transforming the DM from zeeman basis to the eigen-basis

%check diagonalization!!!        
HH;

%DM evolution at polarization field: washing out the coherences during the
%preparation period

      for i=1:64
          for j=1:64
              
              ex=energ(i,i)-energ(j,j);
if ex==0
    rho(i,j)=rho0(i,j);
else
    rho(i,j)=rho0(i,j)*(exp(-2*pi*ai*ex*tp)-1)/(-2*pi*ai*ex*tp)*exp(-2*pi*ai*ex*te);
end              
          end
      end

rhot=rho;

rhot=evec*rhot*evec1;   % now going back to the Zeeman basis; new DM takes into account the evolution at the preparation field


% DM at polarization field is calculated
% now field variation
    Hi = zeros(64,64);


for ifv=1:Nfv-1

    Bint=10000*BotT(ifv,2);
    dt=BotT(ifv+1,1)-BotT(ifv,1);
    
    Hi=Bint*3.0013e2*(ch1*I1z+ch2*I2z+ch3*I3z+ch4*I4z+ch5*I5z+ch6*I6z)/7.04e4;
    Hi=Hi+J12*p12+J14*p14+J32*p32+J34*p34+J52*p52+J54*p54+J26*p26+J46*p46;
       
    
   rhot=expm(-2*pi*ai*Hi*dt)*rhot*expm(2*pi*ai*Hi*dt);  % evolution during switching; Hamiltonian is Hi and time period is dt
    
    
    ifv
end

H0=zeros(64,64);    % Hamiltonian at the detection field; needed to calculate the NMR spectrum

    H0=Bdet*3.0013e2*(ch1*I1z+ch2*I2z+ch3*I3z+ch4*I4z+ch5*I5z+ch6*I6z)/7.04e4;
    H0=H0+J12*p12+J14*p14+J32*p32+J34*p34+J52*p52+J54*p54+J26*p26+J46*p46;
       
% Applying pulse sequence

fi=pi/4;            % flip angle in NMR
    
    fx=zeros(64,64);
    fxx=zeros(64,64);
    
    fx=expm(ai*fi*Iy);
    fxx=expm(-ai*fi*Iy);
    
    [e0,v0]=eig(H0);
    e01=inv(e0);
    
    rhot=e01*rhot*e0;   % DM for PHIP-NMR: converting from zeeman basis to eigen basis
    rhoeb1=e01*Iz*e0;   % DM for equilibrium NMR
    fx=e01*fx*e0;
    fxx=e01*fxx*e0;
    magn=e01*(Ix+ai*Iy)*e0;
    
    AAA=zeros(64,64);
    BBB=zeros(64,64);
    
    % destroying all coeherences in the DM after the spin system arrives to
    % the detection field
                for i=1:64
                for j=1:64
               if (abs(i-j)>0)
                   rhot(i,j)=0;
                   rhoeb1(i,j)=0;
               end
                end
                end
    
ZA=fxx*rhot*fx; % PHIP DM after pulse
ZB=fxx*rhoeb1*fx; % DM for equilibrium after pulse
magnt=transpose(magn);

% Detection: 
                
    for i=1:64
        for j=1:64
AAA(i,j)=AAA(i,j)+real(ZA(i,j)*magn(j,i));
BBB(i,j)=BBB(i,j)+real(ZB(i,j)*magn(j,i));
                end
    end

 
    
    Nomeg=10000;    % number of pints in the NMR spectrum
    omegF=-11;
    omegL=0;
    
    
%     omegF=-7.4;     % staring and final NMR frequency in the spectrum
%     omegL=-0.8;
    
spe=zeros(Nomeg,3);
% now calculating PHIP-NMR (spe(iw,2) and ordinary NMR (spe(iw,3))
    
    for iw=1:Nomeg
    
            w=omegL-(iw-1)*(omegL-omegF)/(Nomeg-1);
            spe(iw,1)=w;
            w1=Bdet*w*3.0013e2/7.04e4;
for i=1:64
for j=1:64
    
    xx=abs(w1-v0(j,j)+v0(i,i));
    k2=1/T2;
    
   spe(iw,2)=spe(iw,2)+AAA(i,j)*k2/(xx*xx+k2*k2);
   
%   real(1.0/(ai*(w1-v0(j,j)+v0(i,i))+1/T2)/sqrt(2*pi));
      spe(iw,3)=spe(iw,3)+BBB(i,j)*k2/(xx*xx+k2*k2);
      
%      /(ai*(w1-v0(j,j)+v0(i,i))+1/T2)/sqrt(2*pi));
end
end

    end

    % plotting the result    
% Plot PHIP NMR
subplot(2,1,1);
plot(spe(:,1),real(spe(:,2)));
title('PHIP');

% Plot Ordinary NMR
subplot(2,1,2);
plot(spe(:,1),real(spe(:,3)));
title('Thermal equilibrium');
    
% plot(spe(:,1),real(spe(:,2)),spe(:,1),real(spe(:,3)));

% saving the result
save propanal-tp_100-te_0.dat spe -ASCII
toc;